── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
✔ dplyr 1.1.1 ✔ readr 2.1.4
✔ forcats 1.0.0 ✔ stringr 1.5.0
✔ ggplot2 3.4.2 ✔ tibble 3.2.1
✔ lubridate 1.9.2 ✔ tidyr 1.3.0
✔ purrr 1.0.1
── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
✖ dplyr::filter() masks stats::filter()
✖ dplyr::lag() masks stats::lag()
ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors
Attaching package: 'kableExtra'
The following object is masked from 'package:dplyr':
group_rows
Attaching package: 'Matrix'
The following objects are masked from 'package:tidyr':
expand, pack, unpack
here() starts at /Users/johngraves/Dropbox/teaching/SMDM-Europe-2023
Loading required package: forecast
Registered S3 method overwritten by 'quantmod':
method from
as.zoo.data.frame zoo
Registered S3 methods overwritten by 'demography':
method from
print.lca e1071
summary.lca e1071
This is demography 2.0.0.9000
Attaching package: 'MASS'
The following object is masked from 'package:dplyr':
select
Attaching package: 'expm'
The following object is masked from 'package:Matrix':
expm
Learning Objectives
Understand options for modeling background mortality based on life tables and parametric mortality models.
Construct a cause-deleted life table to model cause-specific and non-cause-specific death.
Approaches to Modeling Mortality
Standard approach is to draw from life table data:
(We’ll try to replicate these numbers with a discrete time Markov model in a bit.)
Alive-Dead Model
Alive-Dead Model
As a refresher, let’s construct a very simple alive-dead model using life table death probability (qx) values.
Parameterize the model
params =list(t_names =c("lifetable"), # Strategy names. n_treatments =1, # Number of treatmentss_names =c("Alive", "Dead"), # State namesn_states =2, # Number of statesn_cohort =100000, # Cohort sizeinitial_age =0, # Cohort starting age n_cycles =110, # Number of cycles in model. cycle =1, # Cycle lengthlife_table = lt # Processed HMD life table data (2019, USA) )
Parameterize the model
params =list(t_names =c("lifetable"), # Strategy names. n_treatments =1, # Number of treatmentss_names =c("Alive", "Dead"), # State namesn_states =2, # Number of statesn_cohort =100000, # Cohort sizeinitial_age =0, # Cohort starting age n_cycles =110, # Number of cycles in model. cycle =1, # Cycle lengthlife_table = lt # Processed HMD life table data (2019, USA) )
No specific strategies or policies under consideration; we’re just trying to replicate mortality using life table data for a hypothetical cohort of 100,000 infants.
Parameterize the model
params =list(t_names =c("lifetable"), # Strategy names. n_treatments =1, # Number of treatmentss_names =c("Alive", "Dead"), # State namesn_states =2, # Number of statesn_cohort =100000, # Cohort sizeinitial_age =0, # Cohort starting age n_cycles =110, # Number of cycles in model. cycle =1, # Cycle lengthlife_table = lt # Processed HMD life table data (2019, USA) )
Discrete time will be in annual cycles.
Parameterize the model
params =list(t_names =c("lifetable"), # Strategy names. n_treatments =1, # Number of treatmentss_names =c("Alive", "Dead"), # State namesn_states =2, # Number of statesn_cohort =100000, # Cohort sizeinitial_age =0, # Cohort starting age n_cycles =110, # Number of cycles in model. cycle =1, # Cycle lengthlife_table = lt # Processed HMD life table data (2019, USA) )
The lifetable “strategy” will draw on the HMD life table data processed earlier.
Age-Dependent Transition Matrix
We’ll next construct a transition probability matrix using the age-specific probability of death field (qx) from the life table data.
To do so, we need to define a function that, for a given age and cycle length, calculates the probability of death and places this within a transition matrix.
Age-Dependent Transition Matrix
params$mP <-fn_mPt(1:params$n_cycles, params)
Life Table:
lt[1,c("age","qx")]
# A tibble: 1 × 2
age qx
<dbl> <dbl>
1 0 0.00553
Transition Matrix:
Age-Dependent Transition Matrix
Life Table:
lt[1,c("age","qx")]
# A tibble: 1 × 2
age qx
<dbl> <dbl>
1 0 0.00553
Transition Matrix:
params$mP[[1]][,,"lifetable"]
to
from Alive Dead
Alive 0.99447 0.0055322
Dead 0 1
Age-Dependent Transition Matrix
Life Table:
lt[1,c("age","qx")]
# A tibble: 1 × 2
age qx
<dbl> <dbl>
1 0 0.00553
lt[90,c("age","qx")]
# A tibble: 1 × 2
age qx
<dbl> <dbl>
1 89 0.123
Transition Matrix:
params$mP[[1]][,,"lifetable"]
to
from Alive Dead
Alive 0.99447 0.0055322
Dead 0 1
params$mP[[90]][,,"lifetable"]
to
from Alive Dead
Alive 0.87677 0.12323
Dead 0 1
Construct the Markov Trace
sim_cohort <-function(params) { params$t_names %>%map(~({ tr_ <-t(c("alive"= params$n_cohort, "dead"=0)) res <-do.call(rbind,lapply(params$mP, function(tp) { tr_ <<- tr_ %*%matrix(unlist(tp[,,.x]),nrow=params$n_state) })) res <-rbind(c(params$n_cohort,0),res) dimnames(res) <-list(paste0(c(0:params$n_cycles)), params$s_names) res })) %>%set_names(params$t_names)}
Construct the Markov Trace
Output: named list object with full Markov trace for each state.
Construct the Markov Trace
Output: named list object with full Markov trace for each strategy.
\(q[x]/p[x] = A^{[(x+B)^C]} + D \exp[-(E_i \log(x/F_))^2] + G H^x\)
kostaki
Mortality Modeling
Key takaways:
MortalityLaws has a number of mortality models it can fit to life table data, depending on the context and population you are modeling.
You can easily experiment around with different models to see how well they fit your mortality data.
Mortality Modeling
Another nice feature of the package is that each mortality model type has its own defined function to calculate the mortality rate (or probability) given an age and model coefficients.
So rather than assume a constant mortality rate within age bins, you can model the mortality rate directly for a 40.3 year old, or a 90.0001 year old.
Mortality Modeling
Generally speaking, we need three inputs:
age: Ages for lifetable (e.g., 0, 1, 5, …, 100; or 0,1,2,3,…,100).
dx: The number of deaths between exact ages x and x+1.
lx: Number of survivors to exact age x.
Gompertz Model
Because we want to stay general (i.e., model all over the age spectrum), our first attempt will be a Gompertz model.
Gompertz reduces mortality to two parameters: A and B.
Mortalitiy rate at age \(x\):
Gompertz Model
Because we want to stay general (i.e., model all over the age spectrum), our first attempt will be a Gompertz model.
Gompertz reduces mortality to two parameters: A and B.
Mortalitiy rate at age \(x\):
\[
mx = A \exp(Bx)
\]
Gompertz Model
ages <- lt$age[lt$age<100]deaths <- lt$dx[lt$age<100]exposure <- lt$lx[lt$age<100]gom_fit <-MortalityLaw(x = ages, # vector with agesDx = deaths, # vector with death countsEx = exposure, # vector containing exposureslaw ="gompertz",opt.method ="LF2")
Gompertz Model
# Fitted coefficientsgom_fit$coefficients
A B
0.00010772 0.07534542
Gompertz Model
# Fitted coefficientsgom_fit$coefficients
A B
0.00010772 0.07534542
# Hazard rate of death for a 40 year old based on Gompertz model# mx = A exp(Bx)gom_fit$coefficients["A"] *exp(gom_fit$coefficients["B"] *40 )
A
0.0021938
Gompertz Model
# Fitted coefficientsgom_fit$coefficients
A B
0.00010772 0.07534542
# Hazard rate of death for a 40 year old based on Gompertz model# mx = A exp(Bx)gom_fit$coefficients["A"] *exp(gom_fit$coefficients["B"] *40 )
A
0.0021938
# Can simply use the supplied function to calcualte MortalityLaws::gompertz(x =40, par = gom_fit$coefficients)
$hx
[1] 0.0021938
$par
A B
0.00010772 0.07534542
$Sx
[1] 0.97269
Gompertz Model
# Calculate death hazard rate at age 40 from fitted Gompertz model. mx = gom_fit$coefficients["A"] *exp(gom_fit$coefficients["B"] *40 )
Gompertz Model
# Calculate death hazard rate at age 40 from fitted Gompertz model. mx = gom_fit$coefficients["A"] *exp(gom_fit$coefficients["B"] *40 )# Convert mortality rate to probability-log(1-mx)
A
0.0021962
Gompertz Model
# Calculate death hazard rate at age 40 from fitted Gompertz model. mx = gom_fit$coefficients["A"] *exp(gom_fit$coefficients["B"] *40 )# Convert mortality rate to probability-log(1-mx)
A
0.0021962
# Compare with life table probability (qx) at age 40:lt[41,c("age","qx")]
# A tibble: 1 × 2
age qx
<dbl> <dbl>
1 40 0.00208
Gompertz Model
Gompertz doesn’t fit all age ranges well—particularly young ages.
(This is a well-known fact in demography.)
Gompertz Model
plot(gom_fit)
Gompertz Model
If we were to focus only on adults, however, this would be a nice way to go.
While you could draw from any of the aforementioned mortality models, the Heligman-Pollard model tends to fit the entire age distribution reasonably well.
Easy to use: just use HP in lieu of gompertz!
Heligman-Pollard
While you could draw from any of the aforementioned mortality models, the Heligman-Pollard model tends to fit the entire age distribution reasonably well.
Easy to use: just use HP in lieu of gompertz!
ages <- lt$age[lt$age<100]deaths <- lt$dx[lt$age<100]exposure <- lt$lx[lt$age<100]hp_fit <-MortalityLaw(x = ages, # vector with agesDx = deaths, # vector with death countsEx = exposure, # vector containing exposureslaw ="HP",opt.method ="LF2")
Heligman-Pollard
plot(hp_fit)
Heligman-Pollard
Mortality rate for a 40 year old:
HP(x =40, par = hp_fit$coefficients)
$hx
[1] 0.0020409
$par
A B C D E F_
0.000386408 0.021271772 0.107422357 0.001025285 2.871760263 33.201213579
G H
0.000022277 1.102611045
x =40mu1 = A^((x + B)^C) + G * H^x)mu2 = D *exp(-E * (log(x/F_))^2))eta =ifelse(x ==0, mu1, mu1 + mu2)hx = eta/(1+ eta)
Death Hazard at Age 40: Gompertz vs. HP
gompertz(x =40, par = gom_fit$coefficients)$hx
[1] 0.0021938
HP(x =40, par = hp_fit$coefficients)$hx
[1] 0.0020409
Alive-Dead (Modeled Mortality)
Objectives
Fit a parameteric mortality model to life table data.
Use the estimated parameters to sample death times (useful for DES modeling)
Use the estimated parameters as background mortality rates in discrete time Markov.
1. Sampling Death Dates for DES
Difficult to sample death times from life tables due to coarseness of data (e.g., five-year age bins).
Solution:
Fit a parametric mortality model (e.g., Heligman-Pollard)
Construct a survival function from modeled hazards.
Sample a death time from the survival function.
1. Sampling Death Dates for DES
Intuition is easier if we start with a defined quantile (e.g., 25th percentile of survival)
1. Sampling Death Dates for DES
Intuition is easier if we start with a defined quantile (e.g., 25th percentile of survival)
u =0.25
1. Sampling Death Dates for DES
We can feed the quantile (0.25) and coefficients into an inverse quantile function.
Straightforward to do for Gompertz model
u =0.25qgompertz(p = u, shape = gom_fit$coefficients["B"], rate = gom_fit$coefficients["A"])
[1] 70.467
1. Sampling Death Dates for DES
Let’s verify this gets us close to what we see in the life table data.
u =0.25qgompertz(p = u, shape = gom_fit$coefficients["B"], rate = gom_fit$coefficients["A"])
We now have the tools we need to sample death times that closely match underlying mortality data!
For a given simulated patient, sample a uniform (0,1) and feed this number into the inverse quantile function.
1. Sampling Death Dates for DES
We now have the tools we need to sample death times that closely match underlying mortality data!
For a given simulated patient, sample a uniform (0,1) and feed this number into the inverse quantile function.
# Sample death times for 10 individualsu <-runif(10, min =0, max =1)qHP <-approxfun(y =0:115,x =1-exp(-cumsum(HP(x =0:115, par = hp_fit$coefficients)$hx)))qHP(u)
Let’s now incorporate modeled mortality parameters into our Alive-Dead model.
Rather than draw a probability of death from the life table data, we’ll calculate the probability of death at a given age using the model coefficients.
Parameterize the model
params =list(t_names =c("lifetable", # Strategy names"heligman-pollard","gompertz"), n_treatments =3, # Number of treatmentss_names =c("Alive", "Dead"), # State namesn_states =2, # Number of statesn_cohort =100000, # Cohort sizeinitial_age =40, # Cohort starting age n_cycles =110, # Number of cycles in model. cycle =1, # Cycle lengthlife_table = lt, # Processed HMD life table data (2019, USA)gom_fit = gom_fit, # Fitted Gompertz model objecthp_fit = hp_fit # Fitted Heligman-Pollard model object )
Merge these percentages into the underlying life table data and use them to calculate the total number of CVD deaths by age:
lt_ <-# lt %>%mutate(age_ihme =cut(age,unique(c(0,1,seq(0,95,5),105)),right=FALSE)) %>%left_join(ihme_cvd,"age_ihme") %>%mutate(dx_i =round(dx * pct_cvd)) %>%select(age, age_ihme, pct_cvd, dx, # Deaths dx_i, # Cause-specific deaths lx) %>%# Livingmutate(a =ifelse(age==0, 0.152, 0.5)) %>%# time lived by deaths in age groupmutate(age_interval =c(diff(age), NA)) %>%select(age,age_ihme,pct_cvd, lx,dx,dx_i,a,age_interval) %>%# Probability of death in interval = deaths / total living at beginning of interval. mutate(q =replace_na(1-lead(lx) / lx, 1)) %>%# Convert probability to rate. Note that we could also use q/(age_interval - q * (age_interval - a))mutate(m =-log(1-q)/age_interval) lt_ %>%select(age,age_ihme,pct_cvd,dx,dx_i) %>%filter(age %in%c(0,10,25,50,75,98)) %>%kable() %>%kable_styling()
age
age_ihme
pct_cvd
dx
dx_i
0
[0,1)
0.01675
553.221
9
10
[10,15)
0.04440
11.518
1
25
[25,30)
0.05308
103.706
6
50
[50,55)
0.24371
365.721
89
75
[75,80)
0.31064
1997.766
621
98
[95,105)
0.49531
1262.593
625
We next calculate the cause-deleted probability of death (qd) and death rate (md), as well as the cause-specific probability of death (qi) and death rate (mi):
Steps
Cause-specific probability of death (qi): \(q \cdot dx_i / dx\)
Cacluate cause-specific rates by dividing deaths of a given cause into person-years of exposure.
This is equivalent to multiplying the overall rate by the ratio of deaths of a given cause to the total.
Steps
Cause-specific probability of death (qi): \(q \cdot dx_i / dx\)
Cacluate cause-specific rates by dividing deaths of a given cause into person-years of exposure.
This is equivalent to multiplying the overall rate by the ratio of deaths of a given cause to the total.
lt_ <- lt_ %>%# Cause specific probability of death. mutate(qi = q * dx_i / dx) %>%mutate(Rd = (dx - dx_i) / dx) %>%mutate(md = m * Rd) %>%mutate(mi = m - md) %>%mutate(dxd = dx - dx_i)